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ABSTRACT 

A  systematic  method  is  presented  for  deriving  the  Thomas-Fermi 
equation  for  an  atom  and  the  quantum  corrections  from  the  many  "body 
description.  The  novel  feature  of  the  method  is  that  it  does  not 
require  any  a  priori  assumptions  about  the  distribution  of  the 
electrons  in  phase  space  but  shows  instead  that  the  distribution 
which  is  usually  assumed  is  a  direct  consequence  of  specifying  that 
the  many  particle  system  is  in  its  ground  state.   The  procedure 
used  in  the  derivation  is  the  expansion  of  the  mixed  position-momentum 
representation  of  the  Green's  function  in  a  series  of  powers  of -n. 
The  lowest  order  term  is  found  to  correspond  with  the  Thomas-Fermi 
density.   The  form  of  the  higher  order  terms ,   which  are  to  be  considered 
as  corrections  to  zeroth  order  term,  depends  on  the  approximations  made 
in  the  many  body  equations  for  obtaining  the  Green's  function.   This 
paper  deals  only  with  the  Hartree-Fock.  approximation t   but  the  methods 
presented  here  allow  generalization  to  other  approximations  which  can 
include  correlation  effects. 


I.  INTRODUCTION 

The  Thomas-Fermi  model  of  the  atom  developed  historically 
from  reasonable  physical  assumptions  about  the  nature  of  an  atomic 
system.   Because,  however,  the  development  was  not  a  systematic 
derivation  from  the  many  body  Schrodinger  equation,  it  was  difficult 

to  incorporate  corrections  for  certain  higher  order  effects.  Thus, 

2 
the  reasonable  inhomogeneity  correction  introduced  by  von  Weizsacker 

led  to  results  which  were  less  in  accord  with  experimental  results 

than  those  of  the  uncorrected  method.   A  similar  fate  awaited  the 

introduction  of  the  effects  of  the  Pauli  principle  by  Dirac. 

The  success  of  the  Thomas -Fermi  model  however  would  indicate 

that  it  is  soundly  based.   More  recently,  some  work  has  appeared  by 

h  5 

Kompaneets  and  Pavlovskii  and  by  Kirzhnitz  which  would  indicate 

that  the  Thomas-Fermi  model  is  somehow  related  to  the  Hartree  self 

consistent  field  equations  for  the  atom.  More  precisely,  they  have 

shown  that  the  assumption  of  a  local  Fermi  density  is  consistent  with 

a  first  approximation  of  the  Hartree-Fock  equations. 

In  this  paper,  we  give  a  systematic  derivation  of  the  Thomas-Fermi 
model  from  the  many  body  equations  of  an  atomic  system.  We  use  the 

Green's  function  formulation  introduced  by  Schwinger,   discussed  by 

7 
Galitskii  and  Migdal  for  the  quantum  mechanical  many  particle  system 

Q 

in  the  ground  state,  and  by  Martin  and  Schwinger  for  the  quantum 
mechanical  many  particle  system  at  non  zero  temperature.   This 
formalism  for  treating  the  N  particle  system  is  a  time  dependent 
description  which  resembles  in  some  ways  the  time  independent 


description  employing  the  set  of  1,  2,  ...   N  particle  density- 
matrices.   In  it >  the  n  particle  Green's  function 

G  (rnt  ...  r  t  i  r't*  ...  rV)  (1-1) 

is  analogous  to  the  n  particle  density  matrix  and  can  supply  all 

9 
of  the  information  about  the  system  contained  in  the  latter.   As 

is  the  case  for  the  density  matrices,  the  Green's  functions  satisfy 

a  hierarchy  of  equations  in  which  the  n  particle  Green's  function 

is  coupled  to  the  (n-l)  particle  and  the  (n+l)  particle  Green's 

function. 

Explicit  solutions  to  this  set  of  equations  coupling  the 

Q 

various  Green's  functions  are  not  known.   Martin  and  Schwinger  discuss 
a  systematic  approximation  scheme  for  truncating  this  set  of  equations, 
retaining  the  first  n  equations  and  approximating  G  ,  in  terms  of  the 
Green's  functions  for  fewer  particles.   The  simplest  of  these  approxi- 
mations is  the  Hartree-Fock  approximation,  which  results  from 
approximating  the  two  particle  Green's  function  by  the  antisymmetric 
product  of  one  particle  Green's  functions.  We  shall  show  that  the 
first  term  in  the  expansion  of  the  solution  of  the  Hartree-Fock  equation 
in  powers  of  •#  yields  the  Thomas-Fermi  model.   The  Dirac  and 
inhomogeneity  effects  appear  together  in  the  equation  determining  the 
first  non  vanishing  correction  to  the  Thomas-Fermi  model.   Higher  order 

corrections  can  be  generated  in  a  manner  similar,  but  not  identical  to, 

k 
the  method  of  Kompaneets  and  Pavlovskii.   The  fact  that  the  local  densit; 

is  the  density  of  a  Fermi  gas  is  a  consequence  of  the  theory  and  not  an 

additional  assumption  required  to  solve  the  problem. 


The  technique  of  deriving  the  Thomas -Fermi  model  from  the 
many  body  equation  is  useful  "because  it  would  enable  one  to  incorporate 
correlation  corrections  into  the  Thomas-Fermi  model.  We  shall  treat 

this  problem  in  a  forthcoming  paper  by  replacing  the  Hartree-Fock 

,   .  .    10 
approximation  with  the  Gell-Mann  -  Brueckner  approximation. 

In  Sec.  II,  we  state  without  proof  some  useful  properties  of 
the  Green's  functions.   We  also  state  the  first  hierarchy  equation 
satisfied  by  the  single  particle  Green's  function  and  the  Hartree-Fock 
approximation  to  it.  We  transform  this  equation  to  a  mixed  position-momentum 
representation.   This  transformation  introduces  an  infinite  order 
differential  operator  9,  which  acts  on  the  mixed  representation  single 
particle  Green's  function  to  produce  an  equation  completely  equivalent 
to  the  original  Hartree-Fock  equation.   In  Sec.  Ill,  we  make  use  of 
the  structure  of  9  to  expand  the  solution  in  a  series  involving  powers 
of  -fi.  We  then  show  that  the  lowest  order  term  yields  the  Thomas-Fermi 
density.   In  Sec.  IV,  we  discuss  higher  order  corrections.   Sec.  V 
contains  the  proof  of  the  assertions  made  in  Sec.  II.  Most  important, 
however,  it  contains  the  discussion  of  the  analytic  properties  of  the 
exact  Green's  function  which  indicates  clearly  the  connection  between 
the  appearance  of  the  Thomas-Fermi  density  and  the  specification  that 
the  system  be  in  its  ground  state.  The  derivation  of  the  operator  9 
is  to  be  found  in  the  appendix. 

II.  TRANSFORMATION  TO  THE  MIXED  REPRESENTATION 
We  consider  an  atom  or  an  ion  to  consist  of  a  nucleus  with 
charge  Ze  and  infinite  mass  located  at  the  origin,  surrounded  by  N 


electrons,  each  of  mass  m.  When  we  neglect  spin  dependent  forces 
and  all  relativistic  effects  in  our  system,  the  Hamiltonian  is 

H  =  V.H  (r.)  +  i^~.  ,  .  v(r.  -  r.)    .  (2.1a) 

H  (r.)  =  (-n2/2m)V.2  -  Ze2/|r.|     .  (2.1b) 

v(r.  -  r.)  =  e2/|r.  -r.l       .  (2.1c) 


We  have  mentioned  the  n  particle  Green's  function  in  the 
previous  section:  For  the  purposes  of  this  section,  we  shall  need 
only  a  few  of  its  simplest  properties.   These  properties  follow 
directly  from  the  basic  definition  (see  Sec.  V)  and  will  be  stated 
in  the  following  paragraph.   First,  however,  we  must  comment  on 
notation.   If  the  n  particle  Green's  function  is  denoted  by  (l.l), 
then  each  coordinate  r.  and  r.1  specifies  both  a  spatial  location  and 
a  spin  direction.   Consequently,  equality  of  two  coordinates  indicates 
that  the  spatial  locations  are  the  same  and  that  the  associated  spin 
indices  are  identical.   Integration  over  a  coordinate  will  be  understood 
to  imply  summation  over  the  corresponding  spin  index,  and  a  Dirac  delta 
function  of  two  coordinates  will  be  understood  to  be  multiplied  by  the 
Kronecker  delta  of  the  corresponding  spin  indices. 

The  n  particle  Green's  function  contains  all  the  information 

necessary  to  compute  the  expectation  value  of  any  n  particle  property 

Q 
of  the  system.   In  particular,  the  density  of  particles  with  a  given 

direction  of  the  spin  is  given  by 

n(r,t)  =  (-i)G1(rt;rt+)    ,  (2.2) 


where  t  is  a  time  infinite simally  later  than  t.   The  energy  of  the 
system  is  the  expectation  value  of  the  Hamiltonian  and  is  computed 
from  the  one  and  two  particle  Green's  functions  by  the  prescription 

E  =  (-1)  /  dr  Lim  ,    H  (r)  G,  (rt:  r't) 

+  (-i)2  i   /  drndr  v(r  -rj  G0(r.t  rt:r.tVt+).  (2.3) 
»     2  /   *wl  '^   w"i  ^2'   2  ""l  *>*2  **i  "^2 

For  the  system  governed  by  the  Hamiltonian  (2.1),  the  Green's  functions 
satisfy  a  set  of  integro-differential  equations  in  which  the  n  particle 
Green's  function  is  coupled  to  the  (n-1)  particle  and  to  the  (n+l) 
particle  Green's  functions.   The  first  of  these  equations  is 

[  ifi(a/at1)  -  Ho(Xl)  ]  G1(r1tL;  r^) 

+  ij  ^2v(rx-r2)  G^t^r^j  r^ r^) 

=  *  S^-r^)  5(tx-  t[)      .  (2:k) 

The  Hartree-Fock  approximation  results  from  approximating  Q     in  the 
equation  above  by  the  anti symmetrized  product  of  one  particle  Green's 

functions  as  follows: 

G2(1.,2>1,2*)  =   G1(l>l')G1(252')  -  G^l^G^l')  .    (2.5) 

The  functions  appearing  in  (2.k)    and  (2.5)  have  the  spin  indices 
implicit  in  the  coordinates  and  a  spin  summation  implied  by  the 
integration.  When  we  neglect  spin- dependent  forces  of  our  system 
and  choose  a  Hamiltonian  (2.1)  ,   a  single  particle  propagates  without 
change  of  spin.   That  is,  the  single  particle  Green's  function  vanishes 


unless  the  spin  Indices  for  the  ^article   r.  and  the  spin  coordinate 
for  the  particle  r'  are  the  same. 

Substitute  (2.5)  into  (2.4)  and  make  the  spin  indices  and 
summation  explicit.   If  we  use  the  fact  that  the  Green's  function 
for  a  spin  independent  system  is  zero  if  the  two  spin  indices  are  not 
alike  and  is  otherwise  independent  of  spin,  then  we  find  the  only- 
effect  to  be  a  factor  of  two  multiplying  the  term  arising  from  the 
first  product  on  the  right  of  (2.5).  Accordingly,  we  shall  include 
this  factor  of  two,  and  in  the  future  understand  that  there  are  no 
longer  any  implicit  spin  indices  or  summations  in  the  equation  which 
results: 

[ifi(a/&t1)-Ho(5:1)  ]  G^t^'tp 

+  i  J  ^2  v^l-^2)  [2G1^2t'£2t+)  ^(£1^1) 

«*  BCr^-r^)  5(t1-t;J)    .  (2.6) 

It  is  in  the  spirit  of  the  Hartree-Fock  self  consistent  method 
to  consider  (2.6)  as  a  pair  of  integral  equations  by  writing  it  as 

J   K(r1t1jr2t2)dr2dt2G1(r2t2jr^)  =  n  sQ^-r^)  5(^-1^),  (2.7] 
K(£1t1jr2t2)=  [in^/ot^  -  Ho(Xl)  +  2i  J   drVCr-L-  r  '^(r  'tjr  't+)] 
X8(tx-t2)  5^-rg)-  iv(r1-r2)G1(r1t1;r2t1+)5(t1-t^)  .  (2.8; 


Taking  the  Fourier  transform  of  (2.7)  with  respect  to  rn  -  r'  and  t, -  t ' , 

'  ^1  «v>l      1   1' 

keeping  R^  =  p(&,  +  r  )  fixed,  we  have 

/  *(Xl\ ' r2t2)  ^Vltes^  **$)  d(%  '  ^1}  d( *1 '  tl) 

x  exp  j  -i[p  •  (^-r^-^-t^)]  /*>-*.  (2.9) 

The  integral  in  (2.9)  may  be  expressed  in  terms  of  the  Fourier  transforms 
K  and  G, 

K(R,pco)    =    /   d(r,   -r_)    d(tn  -tjK(rnt    -rt   ) 
\w\Sv  '       J      wl    *~2'      v   1       2       "*1  l'~2  2' 

x  exp  j    -±[p,-    CSa.\E2)-«(t1-t2)]/*    r    >  (2.10a) 

5(R,pu>)   -/dC^-JEg)!^-  t2)G1(r1t1;r2t2) 

I  -*[*  '  QEi-JEaJ-^i-V]/*  j  '         (2'10b) 


X  exp 


R  =  J(rn  +r)    ,  (2.10c) 

ynt  2     ~"1        "~2 

by  means  of  a  differential  operator  of  infinite  order  9  ,        as 

9   [k(R,pu),  G(R,pu)l 

x\p  l  w  w       »w  vw  j 

Equation  (2.9)  then  takes  the  form 

G  [k,5]   =-6  .  (2.11) 

K(R,pu)  may  be  expressed  in  terms  of  G  by  using  the  definitions  (2.8) 

VW  v.* 

and  (2.10) ,   with  the  result  that 


R(R,pu>)  =  oJ-p2/2m+  Ze2/|R|  -  2   dr 'v(R  -r ')n(rj) 


+  /'  -^C  n^p'Rp-pJ)  ,      (2.12) 
J   (2,tft)5    """""    ~  ^ 


where 


+      r   dp 

n(R)  =  -iG(Rt:Rt  )  =  /  £-3  n(R,p)  (2.13b) 

~  "*     J    (2rtR)5   ~  ~ 

v(p)  -JV^e"^*^  (2.3A) 


The  relationship  between  n(R,p)  and  n(R)  jjermits  us  to  interpret 
n(R,p)  as  the  phase  space  density  for  electrons  of  a  given  spin 
direction.   The  system  of  equations  (2.11)  and  (2.12)  is  completely 
equivalent  to  the  original  Hartree-Fock  equation  (2.6).  We  next 
turn  our  attention  to  the  expansion  of  the  set  (2.11)  and  (2.12)  in 
a  power  series  in  ft. 

III.  EXPANSION  IN  POWERS  OF  ft 
We  refer  to  the  appendix  and  find  that  the  operator  9   consists 
of  an  infinite  series  of  terms,  each  term  containing  explicitly  some 
power  of  ft  as  a  factor 


The  zeroth  order  term  is  simply  the  product  of  the  transforms  K  and  G 
on  which  Q       operates.  Higher  order  terms  all  contain  derivatives 
of  the  transform  with  respect  to  components  of  R.  Hence,  in  an 


infinite  homogeneous  system  for  which  neither  K  nor  G  can  depend  on  R, 
the  zeroth  order  approximation  to  (2.11)  is  exact,  within  the 
Hartree-Fock  approximation.  For  a  non  infinite  system,  higher 
order  terms  of  ©  do  contribute,  and  provide  corrections  for  the  effects 
of  inhomogeneities.   The  higher  order  terms  of  the  operator  ©  all 
contain  the  operator  -n(d/dR. )  (d/^P.. )  >   which  suggests  that  the 
physically  significant  parameter  exhibited  by  ©  is  the  fractional 
variation  of  K  and  of  3  within  a  cell  of  area  -ft  in  the  R.p.  phase 
plane.  If  "h  were  small  enough  relative  to  the  extent  of  the  phase 
plane  occupied  by  the  system,  then  higher  order  terms  of  the  operator 
©  would  again  be  unimportant. 

This  observation,  that  the  operator  ©  involves  powers  of  -ft 
in  a  series  of  terms  whose  relative  size  depends  on  the  importance  of 
the  uncertainty  principle,  suggests  that  it  might  be  useful  to  expand 
the  density  n(R)  in  a  series  in  which  the  size  of  ft"  enters  in  a 
systematic  manner.   Accordingly,  we  write 

In  writing  this  expansion  we  do  not  imply  that  n . (r)  is  independent 
of  ft;  if  it  were,  this  would  require  that  there  be  a  classical  charge 
density  for  the  atomic  system,  namely,  n  (r) ,  which  could  be  found 
in  the  limit  ft"  =  0.  As  we  shall  see,  the  quantum  effect  which  prevents 
the  atom  from  collapsing  into  the  nucleus  will  be  manifest  in  a 
negative  power  of  ft.  multiplying  each  of  the  terms  n . (R) . 
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The  expansion  of  the  density  n(R)  implies  the  possibility  of 
similar  expansions  for  n(R,p)  and  G-(R,pw)  : 

n(R,p)  =  n^n^n  (R,p)   ,  (3-3) 

G(R,pw)  =  h  7~  -fi'5G.(R,pu)  •  (3.^) 

The  powers  of  ft"  appearing  outside  of  the  summations  ahove  have  been 
so  chosen  that  those  relationships  between  the  coefficients  which 
are  a  consequence  of  (2.13)  will  be  independent  of  n.   These 
relationships  are 

n.(R,P)  =  -i  f   fj  e^^GjCR^uO  (3-5a) 

P  dp 
n.(R)   =   /  _n_  n  (R,p)   .  (3.5b) 

We  should  now  consider  the  effects  of  these  expansions  on  the  kernel 
K(R,pu) .   The  exchange  potential  (2.1*J-)  for  the  Coulomb  force  (2.1c) 

VW  v»J 

is 

v(p_p«)  =  ^e%2/  Ip-p'I2  •  (3-6) 

Substituting  the  expansions  (3.2)  and  (3«3)  into  the  definition 
(2.12)  of  the  kernel  K  gives 

K(R,pw)   =  u  -p2/2m  +  Ze2/R    -  Y~  hJ   2    /   dr  v(R-r)n.(r) 

+     E  4(J«)  Jssl    f  ^Pl^  n  (k,e-)  .  (3.7) 

j=o  (at)9    J     |p-p  I2 


n 


This  form  suggests  that  we  can  classify  the  terms  in  K  according  to 
the  explicit  power  of  ft  they  contain.   Therefore,  we  write 

K(R,pw)  =J~,   tfJK.(R,pu>)  (3.8a) 

K    =  u- jo2/2m+Ze2/|R|  -2    /  drv(E-r)   nQ(r)  (3.8b) 


K-,    =   -2       dr  v(R  -r)    n.  (r)  (3-8c) 

r                                           Uite2        f           dP 
K    =  "2   /   *r  v(R-r)    n  (r)  +    —3     / n       (R,p').    (3.8d) 

j  j  j  (2jtr  1         «i2      J-2      ~ 

p  -p 


«v 


The  expansions  (3.1),  (3.^)  and  (3«8a)  are  to  be  inserted  in  Eq. 
(3.11)  and  the  explicit  powers  of  ft  are  to  be  separated  in  such  a 
way  that  the  zeroth  order  approximation  is  non  vanishing.   The  result 
is  the  series  of  equations  below: 


[k  ,£4g1  =  tf  (3.9a) 

L  o    oJ 


o 


j+k+£=n   J 

The  solution  to  the  set  of  Eqs.  (3«8)  and  (3-9)  can  be 
obtained  in  principle  by  starting  with  the  zeroth  order  equation  of 
each  set,  solving  simultaneously,  then  proceeding  to  the  first  order 
equation  in  each  set,  solving  these  simultaneously  using  the  results 
of  the  zeroth  order  solutions,  and  so  on.  We  shall  now  show  that  the 
solution  of  the  zeroth  order  equations,  i.e.  of  (3.8b)  and(3»9a), 
yields  the  Thomas-Fermi  model. 

The  operator  Q     appearing  in  (3>9a)  consists  of  the  product 


11a 


(x)-jJL 


PATH"  OF  INTEGRATION 


IN  COMPLEX    oj  PLANE 
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of  the  transforms,  so  that  this  equation  may  he  written 

G0(£>^   =  ^"5  [^(S^)]"1  (5.10a) 

=  h~5  [u  -  e]"1  (5.10b) 

where       E(R,p)  =  P  /2m  +  0  (R)  (5-lOc) 

0o(R)  =  -Ze2/|R|  +  2  J  v(R-r)  nQ(r)dr  .       (5.11) 

It  is  necessary  to  recover  n  from  G  by  performing  the  integrations 

o      o 

indicated  by  Eqs.  (5-5).   In  order  to  carry  these  out,  however,  it 
is  first  necessary  to  specify  a  path  in  the  complex  w  plane.   This 
choice  of  path  is  discussed  in  Sec.  V.   It  turns  out  that  the  choice 
of  a  path  is  determined  by  the  energy  state  of  the  ion  or  atom,  and 
by  the  number  of  electrons  N  surrounding  the  nucleus.   For  a  system 
in  its  ground  state,  the  path  of  integration  lies  just  below  the  real 
w  axis  from  -co  <  w  <  u,  crosses  the  real  axis  at  oj  =\i   ,  and  lies  just 
above  the  real  axis  from  u  <  w  <  +oo  .   The  constant  \i   is  the  chemical 
potential,  approximately  the  energy  required  to  add  another  electron 
to  the  system.  It  determines  the  number  of  electrons  surrounding  the 
nucleus.   This  point  is  brought  out  in  detail  in  Sec.  V. 

The  exponential  e     in  (5- 5a)  makes  it  possible  to  complete 
the  path  by  adding  the  integral  over  the  infinite  semicircle  In  the 
upper  half  of  complex  w  plane,   n  (R>p)  is  then  given  by  a  contour 
integral,  and  is  equal  to  -n   or  zero,  depending  on  whether  the  point 


U) 


=  E  is  inside  or  outside  the  contour.   Since  the  contour  crosses  the 
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real  axis  at  u>  =  u. ,   this  means  that 

nQ(R,P)  =  *"3  j     p2/2m  +  0Q(R)  <  n  (5.12a) 

=  0        p2/2m  +  0  (R)  >  u  .  (3-12h) 

The  density  n  (R)  defined  by  (5«5h)  is  then 

n  (R)  =  (eAVS^R)  (3.13) 

O  r 

where  p  (R)  ,  the  maximum  momentum  of  electrons  in  the  Fermi  gas,  is 
F 

given  by 

PF2/2m+  0Q(R)  =  u  .  (3.1*0 

_■* 

The  density  n  (R,p) ,   which  equals  -ft  J   for  p  <  p_(R) ,   is  of  course 
o  a 

the  density  usually  assumed  at  the  outset  in  derivations  of  the  Thomas- 
Fermi  model.  Its  appearance  here  is  a  consequence  of  the  theory. 

To  complete  the  derivation  of  the  model,  one  may  obtain  the 
integral  form  of  the  Thomas-Fermi  equation  by  eliminating  0  (R)  from 
(3.14)  and  (3.11),  and  eliminating  n  (R)  from  (3«H)  with  the  result 

PF2(R)/2m  -  Ze2/|R|  +  (jA5)'1/  Pp5  (r)  v(R-  r)dr  =  u  .   (3.15) 

The  differential  form  of  the  Thomas-Fermi  equation  may  be  obtained  by 
first  combining  (3.13)  and  (5.1*0  as 

3/2 


n  (R)  =  (6A5)"1 


2m 


[n-0o(R)]  \  '    '  ,  (3-16) 


then  taking  the  Laplacian  of  (5.H)  with  v(R-r)  =  e  /  I R  -  r  I  -   The 


Ik 


result  of  this  operation, 

V20q  =  -8rte2n0(R)   , 
is  combined  with  (3«l6)  to  give 

V20q(R)  =%(^e2/5rtfi5)   2m  [u-  0Q(R)]      ,   (3.17) 
which  is  the  Thomas-Fermi  equation  in  its  usual  form. 

IV.   HIGHER  ORDER  TERMS 
Having  found  that  the  zeroth  order  term  in  the  spatial  density- 
is  the  Thomas -Fermi  density,  we  turn  our  attention  to  a  consideration 
of  the  higher  order  terms.  We  shall  find  that  the  first  order  correction 
vanishes,  and  that  the  second  order  correction  is  identical  with  that 
discussed  "by  Kompaneets  and  Pavlovskii.   We  consider  the  n=l  equation 
of  the  set  (3.9b).   The  operator  9,  (see  appendix)  constructs  the 
Poisson  "bracket  of  the  functions  on  which  it  operates  -  in  this  case 

K  and  G  .   The  zeroth  order  result  (3.10a)  guarantees  that  this 
o      o 

Poisson  bracket  vanishes.  Moreover,  since  the  operator  0  merely 
forms  the  product  of  the  functions  on  which  it  operates,  the  solution 
of  the  first  equation  of  (3-9"b)  is 

Gl  =  -KlGo/Ko  =  -*\/^  "E)2  •  ^•1) 

The  density  n  (r)  may  be  recovered  from  G..  by  the  integration 

indicated  in  (3«5) 

+ 

w    °       ( u  -E) 
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We  have  previously  found  that,  over  the  path  described  in  Sec.  Ill, 

+ 

/,    XUO 


S(x)  =1     x  >  0  (^.5b) 

=0     x  <  0  .  (^«5c) 

Repeated  differentiations  of  (h.J>)   with  respect  to  E  will  develop  a 
series  of  relationships  we  shall  use  later  in  this  section. 

/   did  e  l/o 


2jti(u  -E) 


»      ^hy.[h)      8<"-E>         <*•*> 


Furthermore,  if  polar  coordinates  are  introduced  into  momentum  space, 
one  may  write 

dp  =  m  p  dE  d0  dcosQ  (^-5) 

The  integral  in  (4.2)  may  be  evaluated  using  (k*k)    and  (4.5) •  We  note 
that  when  E  =  n,  p  = 

n^R)  =  (2rt%5)-1mpF(R)  K^R)   -  (4.6) 

An  integral  equation  for  K,  (R)  results  by  using  this  density  in  the 
definition  of  K_(R)  given  by  (3.8c) 

Kj_(R)  =  -(AV^/irp^r)  KL(r)v(R-<r)  (4.7) 

The  integration  in  (4.7)  is  over  that  region  for  which  p„(R)  is 
greater  than  zero  because  the  density  n,  (R)  given  by  (4.6)  is 
non  zero  for  positive  p  only.  Depending  on  the  choice  of  u ,  the 
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region  of  integration  may  or  may  not  be  all  space. 

One  solution  to  (^.7)  is  that  K,  (R)  vanishes  identically. 

To  show  that  this  is  the  only  solution,  multiply  (^.7)  by  K1(R)p1?(R) 

and  integrate  over  the  region  for  which  p_(R)  is  greater  than  zero. 

J? 

Making  use  of  the  specific  form  for  v(R  -  r) ,   the  result  is 

where 

C  =  me2/(rt%3) 

Pfe?   =  h®**®    ' 

The  integral  on  the  left  is  non  negative.  The  integral  on  the  right 
is  the  self  energy  of  a  distribution  of  charge  with  density  pj  it  too 
is  non-negative.   The  minus  sign  then  requires  that  each  integral  vanish. 
Equation  (^.1)  indicates  that  G,  also  vanishes. 

We  may  now  consider  the  n=2  equation  of  (3»9h) ,   which  we  solve 
explicitly  for  G  using  the  fact  that  K,  and  G,  are  zero. 

G  =  -K-G  /K  -  9n    Tk  ,G  1/K 
o    ^2.  o'    o         2  Lo'  oj'  o 

or,  equivalently, 

-*\  =  K2(R,p)/(u)-E)2  +  (W-E)_1  92[(W-E),(W-E)"1]  .  (4.8) 

When  the  set  of  differentiations  which  constitute  6  are  performed,  the 
result  is 
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-fi5G0  =  (oj-E)_2Kc(R,p)  +  (w-E)"3C(R)  +  (u  -  E)  ~Vr,p)   0-9a) 


c(R) 

1 
"51 

d&p) 

l 

~  4 

d20    0  d0 


dR 


R  dR 


d  0   /p  •  R\ 


.   d0   /pxR\  2   ,   /d0 


0-9b) 
2-i 


dR 


mR 


R  dR 


mR 


+  m  VdR 


.  (l+.9c 


The  density  n  (r)  is  to  be  recovered  from.  G     by  integration,  as  before. 
In  performing  this  integration,  we  note  first  that  equation  ($.8h) 
indicates  that  K  will  depend  only  on  the  magnitudes  of  R_  and  p. 
Utilizing  (k.k)   and  0.5)  and  integrating  by  parts  where  necessary, 
we  obtain 


k&\{§)   =  amp^R^)  -  mC(R) 


dp 
dE 


E=u 


m 

+  z 


oE^ 


p  /  D(R,p)  d  cos  9 


2n 


jE=u 


O-10) 


The  remaining  integration  in  (^.10)  may  be  carried  out  by  choosing 
the  polar  axis  in  the  direction  of  R,  after  which  the  differentiations 
may  be  performed,  resulting  in 


l^A5n2(R)  =  2mpF(R)  K2(R,pp)  -  2m2c(R)  /(3PF) 


m 


(Sr)2/<V>   • 


(it.ii) 


A  single  integral  equation  for  KU  will  result  if  the  equation  (3»8d) 
for  Kp  is  written  using  the  density  n  given  by  O-ll)  and  n  given  by 
(3.12).  Instead  of  doing  this,  however,  we  shall  exhibit  the  differentia] 
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equation  which,  shows  that  the  second  order  term  in  our  solution  is 

the  same  as  the  first  quantum  correction  found  by  Kompaneets  and  Pavlovskii. 


Define 


02(R)  =  2  J    n2(r)v(R-r)dr 


(4.12) 


Since  v(R  -  r)  is  the  coulomb  potential, 


V20Q  =  -8ite2n, 


(*.«) 


The  nQ  of  Eq.  (3-12)  used  in  the  j=2  term  of  Eq.  (3.8d)  gives 


K^(R,p)  =  -02  +  fz      J 


ft~5dp 


P-P 


«|2 


P  <P, 


■02(B)  +  ^ 


■F 


2n 


PF-P 


PF+P 


.  (luiio 


Finally,  n  is  eliminated  from  (4.11)  using  (4.13)  and  KL(R,pv)  is 

eliminated  from  (4.11)  using  (4.14).   In  the  equation  which  results, 

-1  2 
we  evaluate  C(r)  using  (4.9b)  as  (4m)   V  0  and  p„(R)  using  (3«l4) 

o  F 


as    [2m(n-0o)]1/2 


V%  +  ^    [2m(,  -  0J]1/2  02  «  <k2e\  A6)"V  -  0O) 


me2r2m(u-0  n-1/2      _     Q  .   /d0 N,    2  _ 


This  equation  was  derived  by  Kompaneets   and  Pavlovskii  using  a  different 
method  in  which  the  Thomas-Fermi   distribution   (3-12)    is  assumed  at  the 
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outset.  The  authors  discuss  this  equation  and.  its  numerical  solution 
and  point  out  that  the  exchange  effect  (first  term  on  the  right)  and 
the  inhomogeneity  correction  (last  term  on  the  right)  appear  together 
in  this  equation.  Therefore  the  first  non  vanishing  correction  to 
the  Thomas-Fermi  equation  contains  both  exchange  effects  and 
inhomogeneity  effects;  this  helps  to  explain  why  attempts  to  include 
one  effect  or  the  other  effect  singly  in  the  Thomas-Fermi  model  led 
to  results  so  far  from  accord  with  experiment. 

V.  GENERAL  PROPERTIES  OF  THE  GREEN'S  FUNCTION 
In  this  section,  our  main  concern  is  to  establish  the  connection 
between  the  state  of  the  ion  or  atom  and  the  path  of  integration  used 
to  recover  the  density  from  the  Green's  function.  This  connection 
follows  from  certain  general  properties  of  the  Green's  function.   These 
properties  are  best  exhibited  by  starting  with  the  well  known  definition 
of  the  Green's  function  and  working  in  the  formalism  of  second 
quantization.   In  this  formalism,  the  Hamiltonian  (2.1a)  takes  the  form 


H  =  Z   /  i|r  +(rt)H  (r)i|r  (rt)dr 
a  J   a  ^   o  «•  a  '« 


+ 


^2    /  i   +(rt)t  t  (r't)v(r-  r')t  ,(r't)^  (rt)drdr'  (5.1) 
2  aa'  J   a  ~"   a'  ~~     «•*  —   a'  ~*    a  ~~  '  —  <~  v^   ' 


if   (rt)  and  \|r   (rt)  are  the  annihilation  and  creation  operators  for  the 
electrons;  a  is  the  spin  index  and  takes  on  two  values.   These  operators 
satisfy  anticommutation  relationships 
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*a(£0  ta«(r,'t)  +  ta,(r't)  lra(rt)  =  0 

It  will  be  convenient  to  suppress  the  spin  index  a  by  considering  it 
to  be  implicitly  contained  in  r  and  using  the  convention  described  in 
the  beginning  of  Sec.  II.  We  shall  use  the  Heisenberg  representation, 
in  which  the  time  dependence  of  the  operators  is 

.f.y.  iHt/fi   ,,      .    -v      -iHt/n 


(5-3) 


,  +  /    ..  iHt/n   ,  +  /      .    -x      -iHt/n 

t   (Ew   =   e  ty    (r  t=0)    e  ' 


States  of  the  ion  are  specified  by  N,  the  number  of  electrons 
surrounding  the  nucleus,  the  energy  E,  and  possibly  some  degeneracy 
parameter  7 ,  i.e., 

H|N,E,7>  =  E|NE7>  (5.^) 

NQp  =  |N,E,7>  =  N|N,E,7>  •  (5ATd) 

The  operator  giving  the  number  of  electrons  surrounding  the  nucleus  is 

Kop  =  jV(rt)  t(rt)  dr  (5-5) 

and  hence,  the  operator  n   (rt)  giving  the  spatial  density  of  electrons 
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with  a  given  direction,  of  spin  is 

n   (rt)  =  \|r+(rt)  \|r(rt)  .  (5.6) 

We  shall  he  especially  concerned  with  the  ground  state  of  the  system 
with  N  electrons,  having  energy  E  (N) ,  and  for  this  state,  we  shall 
suppress  the  indices  NE7  and  the  vertical  bar,  leaving  just  the  angular 
bracket  to  denote  the  state. 

Green's  functions  for  the  many  particle  system  are  well  known 

T  R    Q 

and  have  been  widely  studied. ' '  '*     To  fix  the  definition  we  have  been 
using,  the  n  particle  Green's  function  for  the  ground  state  of  the 
K  electron  system  was  defined  by 

G  (1,2, ...n;  l'.-.n1) 


n 


(-i)\(Kl)^(2)...t(n)/(n')...t+(l'j)+>  (5.7) 


The  operators  in  (5-7)  a**e  the  time  dependent  annihilation  and 
creation  operators;  the  bracket  f   J   indicates  the  time  ordering 
of  the  operators  enclosed:  operators  with  later  time  coordinates 
stand  to  the  left  of  operators  with  earlier  time  coordinates  and  the 
collection  is  multiplied  by  +1,  the  parity  of  the  permutation  of  the 
operators  from  the  sequence  indicated  in  (5-7)  to  the  time  ordered 
sequence.   The  expectation  value  is  taken  with  respect  to  the  ground 
state  of  the  N  electron  system. 

It  is  now  evident  that  Eqs.  (2.2)  and  (2.3)  are  nothing  other 
than  the  statements  that  the  energy  of  the  ground  state  is  equal  to 
the  expectation  value  of  the  Hamiltonian  in  the  ground  state  and  that 


22 

7,8,9 


the  particle  density  is  the  expectation  value  of  the  density  operator 
The  equation  (2.k)    linking  the  one  particle  and  two  particle  Green's 

Q    Q 

function  may  be  derived  '  by  differentiating  (5.3) 

i*dt(£t)/dt  =  iKrt)H-HV(rt)  >  (5-8) 

evaluating  the  commutator  on  the  right  of  (5*8)  using  the  relations 
(5-2),  multiplying  the  result  by  -i\|r  (r't')  from  the  right  and  time 
ordering  the  result.  A  delta  function  arises  because 

|:(t(rt)t+(r't<))+=  Q|  f+(r «t»))+  +  S(r-  r')  B(t -t  •)  . 

Finally,  taking  the  ground  state  expectation  value  results  in  Eq..  (2.I4-) 

We  are  especially  concerned  with  the  Fourier  transform  G 
of  the  one  particle  Green's  function  and  its  w  dependence.   The 
transform  was  defined  by  (2.10b),  which  we  rewrite,  using  the  specific 
definition  (5«T)  as 

G(R,a;pu,)  =  -1  j  d(x  -yj d(tx -  ty)<  Va(x,tx)  ^,ty)> 

Xexp|-i  ^(x-yJ-KVty)]^! 

xexp  -1  fc^-^x-y]/*! 


x   y 


+  x 

't  <  t 

x   y 


R  =  i(x  +  y)  .  (5-9) 
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There  is  no  spin  summation  implied  in  (5*9) •   The  time  dependence  of 
the  matrix  elements  is  exhibited  by  using  a  resolution  of  the  identity 
between  the  two  operators  in  each  matrix  element,  and  the  Heisenberg 
time  dependence  (5*3) 


CM.^<V&)|N,^><N,E?'|  C^> 


exp  \    i(tx-ty)[EQ(N)  -E]/h 


(5-10) 


with 


*„(s)  -  +„.(£»*   =  o) 


a. 


+  , 


a  ~^  x 


*  (y)  =  If  (y,t  =  o) 

CG  •*-*-     co  *^  y 


The  summation  in  (5. 10)  is  symbolic,  and  implies  integration  over  the 
continuous  range  of  eigenvalues  also.   The  only  intermediate  states 
which  will  contribute  are  those  for  which  ^N*  |\|/  |n^  does  not  vanish, 
i.e.,  states  where  Nl  =  N+l.   The  energy  E  appearing  here  specifies 
energy  levels  of  the  system  with  (N+l)  particles,  and  will  be  denoted 
E(K+1).   Similarly,  for  the  other  matrix  element  in  (5*9)  >  "the  same 
technique  yields 

^^VV^x^   -  HE;7<ta+(y)|N-l,E,7><N-l,E,7|^(x)> 


X  exp 


i(tx-ty)[E(N-l)-Eo(N)]/h 


.  (5.11) 


The  spatial  Fourier  transforms  are  go  independent,  and  will  be  denoted 


as 


2k 


Hy    I   a(x-y)<*  (x)|m-l,B,7><IH-l,E,7|  ^„+(y)>  exp[~-ip.(x -y)/ft|l 

=  ft"1  A  f"R,ajp,E(N+l)J  (5.12a) 

Z7  /  d(x-y)  <i|r +(y)|N-l,E,7><N-l,E,7|  *  (x)>  expf-ip.(x-y)/«n 


=  ft"  B  JR,ajp,E(N-l)"|  .  (5.12b) 


When  (5.10),  (5«H)  and  (5«12)  are  inserted  into  (5-9),  the  remaining 
time  integration  does  not  converge  and  is  therefore  ambiguous.   The 
integration  can  be  made  unambiguous  in  the  usual  manner  by  inserting 
a  convergence  factor  exp  -T}|t  -t  |   into  the  integrand  and  taking  the 
limit  T]  -» 0  after  the  inversion  to  the  time  dependent  form  is  performed. 
The  result  of  the  integration  with  the  convergence  factor  inserted  is 

a[r,o,;p,E(N+1)1      ^_   B"|R,ajp,E(N-l)l 

G(R,a;pw)  =£ = +}_   ~ = ~ •  (5-13) 

""   ~      E(N+1)  -E  (N)  -  iT) -w       E  (N)  -E(N-l)  +  ±r\  -  w 

The  sums  in  (5*13)  are  over  the  energy  levels  of  the  (N+l)  particle 
system  and  of  the  (N-l)  particle  system. 

At  this  point,  we  call  attention  to  the  fact  that  E  (W)  appears 
in  (5.15)  because  the  state  which  appeared  in  the  definition  of  the 
Green's  function  was  the  ground  state  of  the  system  with  N  particles. 
We  can  make  use  of  this  by  introducing  the  excitation  energy,  a 
positive  quantity  defined  by 

e(N)  =  E(N)  -  E  (N)  >  0  .  (5.1^) 
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It  is  also  convenient  to  introduce  the  chemical  potential  |i,  defined 
as  the  minimum  energy  required  to  add  another  electron  to  the  system, 
i.e. , 

u(N)  =  EQ(N+1)  -  EQ(N)   .  (5-15) 

These  two  definitions  bring  (5. 13)  to  the  form  which  exhibits  its  w 
dependence  most  clearly,  namely 

A[e(N+l)]  B[e(N-l)] 


G  = 


;(N+1)  u(N)  +  e(N+l)  -in  - u    e(W-l)  u(N-l)  -  e(N-l)  +  li\  -  u 


(5.16) 

We  observe  that  the  singularities  of  (5.I6)  in  the  complex  u 
plane  are  poles  and  branch  cuts  adjacent  to  the  real  axis.   The 
singularities  below  the  real  axis  are  contributed  by  the  first  sum 
in  (5.I6)  and  extend  from  u(N)  to  infinity.   The  poles  arise  from 
the  discrete  (bound  state)  energies  of  the  (N+l)  particle  system. 
The  branch  cut  arises  from  the  continuum  state  energies  of  the  (N+l) 
particle  system,  for  which  the  sum  in  (5.I6)  becomes  an  integral. 
The  bound  state  energies  (if  any)  are  expected  to  lie  below  the 
continuum  state  energies  -  i.e.,  have  lesser  excitation  energies. 
Therefore  the  poles,  if  any,  are  expected  to  lie  to  the  left  of  the 
branch  cut.  No  matter  how  many  bound  states  of  the  (N+l)  particle 
system  there  are,  all  the  singularities  below  the  axis  lie  to  the 
right  of  n(N). 

The  singularities  above  the  real  axis  are  contributed  by  the 
second  term  In  (5.16) ,   and  extend  from  u(N-l)  to  minus  infinity.  This 
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time  the  branch  cut  lies  to  the  left  of  the  poles ,  if  any.   All 
singularities  above  the  real  axis  lie  to  the  left  of  u(N-l).   The 
reason  the  two  sets  of  singularities  terminate  at  the  u  is  that  the 
ground  state  of  the  W  particle  system  was  used  to  define  the  Green's 
function.   If  the  Green's  function  had  been  defined  relative  to  an 
excited  state,  then  each  series  of  singularities  would  extend  past 
its  u  a  distance  equal  to  the  excitation  energy  of  this  state.   This 
statement  is  most  easily  verified  by  substituting  for  E  (N)  in  (5. 13) 
the  energy  of  the  state  used  to  define  the  Green's  function. 

In  the  inversion  of  (5.16)  to  recover  the  time  dependent  form, 
the  path  of  integration  of  to  along  the  real  axis  defines  a  path  which 
passes  beneath  one  set  of  singularities  and  above  the  other.   The 
relationship  between  the  two  sets  of  singularities  depends  on  whether 
u(W)  is  larger  or  smaller  than  u(N-l).   The  exact  relationship  between 
these  two  will  depend  in  detail  upon  the  nature  of  the  system.   For 
an  ion,  however,  the  fact  that  it  becomes  increasingly  difficult  to 
remove  successive  electrons  indicates  that  u  increases  with  N.   Hence, 
the  set  of  singularities  below  the  axis  lies  completely  to  the  right 
of  the  set  of  singularities  above  the  axis,  with  no  overlap. 

The  non- overlapping  of  the  two  sets  of  singularities  makes  it 
possible  to  shift  slightly  the  path  of  integration  off  the  real  axis, 
so  that  it  still  passes  beneath  the  one  set  of  singularities  and  above 
the  other,  even  as  the  operation  of  letting  r\   go  to  zero  moves  all  of 
the  singularities  onto  the  real  w  axis.   In  this  way,  the  imaginary  part 
of  the  denominators  in  (5 .16)  may  be  discarded  but  the  path  used  for 
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recovering  the  time  dependent  form  now  lies  below  the  real  axis  from 
minus  infinity  to  u,  crosses  the  real  axis  at  w  =  jj.  and  lies  above 
the  real  axis  from  u.  to  plus  infinity.  Here,  u  is  some  real  number 
such  that 

li(K-l)  <  u  <  n(N)  (5-17) 

The  path  of  integration  discussed  above  is  possible  for  the 
ground  state  Green's  function,  because  of  the  separation  of  the  sets 
of  singularities.   But  this  path  is  precisely  the  one  which  in  Sec.  Ill 
yielded  the  Thomas-Fermi  density. 
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APPENDIX  -   THE   OPERATOR  9n 

R,P 


An  operator  G^   was  defined  in  Sec.  II  by  the  requirement 
-RjP 

that  it  produces  the  Fourier  transform  of  the  matrix  product  of 
two  functions  when  it  operates  on  the  transforms  of  the  individual 
functions.  The  precise  definition  is 

r~      «     -,   r ■  -ifp-CR-y)  -«(t  -t  )"]/« 

^pL^P"),  G(R,P")  J  =J   e   ^  X       yJ   d(x-y)d(tx-ty) 


X  /  K(xt  :wt  )dwdt  G(wt  :yt  )  .         (A.1) 


WW       £1   **V  "W 


A  four  dimensional  notation  will  be  convenient,  in  which 

px  =  p  •  x  -  tot 
dx  =  dxdt 
dp  =  dpdw 

yW 

R  =  ^(x+y)  • 

The  kernels  K  and  G  in  (A.1)  must  be  expressed  in  terms  of  their 
transforms  so  that  the  right  side  of  (A.l)  becomes 

a    f  -ilp(x-y)  -p  (x-w)  -p p(w-y)|/« 

(2rth)"°J  d(x-y)dw  dpx  dp2  e   L 

,x+w   .  ~  ,w+y 


XK  (^,PX)  3  (^L,p2)   •  (A. 2) 


29 


Eliminate  p  and  p£  by  putting  p  =  p  +  q   p  =  p  +  q^. 
The  preceding  expression  (A. 2)  becomes 

o    r  iflnO- w)  +qjw-y)l/n 

(2rth)        /   d(x-y)dw  dq^q^e    L  x  J 

XK(R+   ^jp  +  q1)    SCR-^jp+qg)      .  (A. 3) 

Each  of  the  transforms  K  and  G  may  be  expanded  as  a  Taylor's  series 

about  R  and  p.  For  clarity  we  shall  temporarily  neglect  the  fact  that 

R  and  p  each  have  four  components.   Then  the  Taylor's  series  expansion 
is 


The  powers  of  (w-y)and  of  (x -w) are  derivatives  of  the  exponential 
with  respect  to  iq, /n  and  iq^/fi.   On  integration  by  parts,  the 
derivatives  of  the  exponential  may  be  transferred  to  the  other  factors 


of  the  integrand,  which,  exclusive  of  the  exponential,  is  then 


jklm 


The  integration  may  be  performed  first  over  w,  then  over  q^,  then 
over  (x  -y)  and  finally  q  ,  with  the  result  that  the  integral  equals 
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the  value  of  (A.k)   with  q.  =  0,  q_  =  0.   Setting  the  q.  equal  to  zero 
deletes  all  terms  from  the  summation  except  those  for  which  i  =  k  and 
m  =  J .  The  result  is 

(a.5) 

The  modification  of  (A.5)  brought  about  because  of  the  multicomponent 
nature  of  R  and  p  is  simply  the  appearance  of  an  index  and  numerical 
factor  for  each  component  -  i.e., 


q^p[k(r,p),g(r,p)] 


r<*  /0xni+Vn3+ VWVH 

-  V"  V  ,\lfl/5J 

2_  2_  n .  in.  in7  Jn.  »k    !k-  »k,  !k,  ! 

nln2n3ni+  klk2k3ki+  1^3^123^ 


-".^^•©'•fer-fef'PiP-ftT'5 


(A.6) 


R,    =  -(t    + 1   )  p,    =   -co 

k       2>   x       y'  ^4 


There  will  be  no  R,  dependence  in  either  of  the  functions  K  or  G, 
hence  only  terms  with  n,  and  k,  both  equal  to  zero  will  contribute. 

If  we  separate  out  terms  according  to  the  explicit  power  of  -n  as 

12 
indicated  by  (3«l)>  then 


oM  - 


KG 


1   .    v-  dK  dG         dK     dG 


1  .  s2 


g2[k,g]  -     (^i) 


xyz 


H 


dx  dp  dp.,  ^x 


:•: 


1 


'b2*.     b 


-g      d2^  s"g 


2  U  2  n  2    .  2  -,  2 
xyz   \dx   dpx    dpx  dx 
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(A.T) 
(A.8) 


y  (  d^K   d2G   +  d^   d2G 

^z  I  s^F  ^pxspy     SpxdPy  ^d7 


*—        *—      dxdp    "dydp" 


P^yP,     xyz 


■y 


(A.9) 
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